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SUMMARY 

The flux-vector and flux-difference splittings of Steger- Warming, Van Leer and Roe 
are tested in all possible combinations in the implicit and explicit operators that can 
be distinguished in implicit relaxation methods for the steady Euler and Navier-Stokes 
equations. The tests include one-dimensional inviscid nozzle flow, and two-dimensional 
inviscid and viscous shock reflection. Roe’s splitting, as anticipated, is found to uniformly 
yield the most accurate results. On the other hand, an approximate Roe splitting of the 
implicit operator (the complete Roe splitting is too complicated for practical use) proves 
to be the least robust with regard to convergence to the steady state. In this respect, the 
Steger- Warming splitting is the most robust: it leads to convergence when combined with 
any of the splittings in the explicit operator, although not necessarily in the most efficient 
way. 


1. INTRODUCTION 

Over the past five years conservative upwind differencing methods have been yielding 
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impressive results for compressible aerodynamics problems[l-4]. Among the techniques 
available for achieving the upwind bias the most prominent and most used ones are the 
flux-vector splittings (FVS) by Steger and Warming[5] and Van Leer[6], and the flux- 
difference splittings (FDS) by Roe[7] and Osher[8]. 

Until recently, most research efforts have been concentrated on implementing the dif- 
ferent splitting formulas in the explicit operator by which the spatial residual is computed. 
Some comparisons have been made[2,9,10], from which the following conclusions can be 
drawn. 

(i) Van Leer’s flux-vector splitting by design and in practice is more accurate than Steger- 
Warming’s, for both smooth and shocked flow[2]. 

(ii) Flux-vector splitting is less accurate than flux-difference splitting because it ignores 
contact discontinuities and shear waves. This inaccuracy particularly shows up as 
excessive numerical diffusion across a boundary layer[9,10]. 

Direct comparisons between Roe’s and Osher’s splittings are not known, but the methods 
should be very close in accuracy. There is a preference for the Roe splitting because of its 
algebraic simplicity, which makes its extension to real-gas physics straightforward[ll], and 
its computational efficiency. For an ideal gas though, the Osher splitting can be coded 
almost as economically! 12]. 

Little research has been devoted so far to the analysis and comparison of these split- 
tings when used in the implicit iterative operators of implicit marching schemes. Particu- 
larly urgent are questions regarding their stability and convergence properties. Barth] 12] 
performed a fixed-point analysis of some linearization of a fully implicit method for the 
one-dimensional Euler equations, incorporating either the exact or an approximate Jaco- 
bian matrix of Roe’s splitting formula(hereafter referred to as RS and ARS respectively). 
Liou[14] made an eigenvalue analysis of the system resulting from the use of the true Jaco- 
bians of the Steger- Warming and Van Leer splittings (hereafter referred to as SWS and VLS 
respectively) in various relaxation procedures for the two-dimensional Euler equations. 

The present paper addresses the following questions: 

Is there an optimal combination of implicit and explicit upwind 
operators, simultaneously achieving stability, accuracy and efficiency 
for a wide class of problems? 

The choice of splitting has been limited to the Steger- Warming, Van Leer and Roe splitting 
formulas. There is no need to use the same splitting in both implicit and explicit operators, 
so there are nine combinations to choose from. Clearly there is a need for a shopping guide. 
The following matters need special attention in this investigation. 

(1) Uniqueness: it is not clear whether the same solution is obtained with different split 
fluxes because of their nonlinear nature. Specifically, one might ask if different implicit 
operators will drive toward the same solution if the same explicit operator is used. 
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(2) Sensitivity: the performance of the implicit operator may be sensitive to the changes 
of residual differencing accuracy, in particular from O(Ax) to 0(Ax 2 ). 

(3) Monotonicity: modern, monotonicity preserving(non-oscillatory) interpolation tech- 
niques for obtaining interface fluxes are known to affect convergence. 

(4) Diffusivity: all splitting formulas ultimately are intended for the computation of the 
inviscid fluxes in the Navier-Stokes equations. The level of artificial diffusion in the 
explicit operator to a great extent determines the accuracy. 

(5) Grid nonuniformity: grid stretching tends to slow down convergence; how this effect 
varies among the different flux formulas is not known. 

The emphasis will be on (l),(2),(3) and (4). Items (4) and (5) often appear together 
in Navier-Stokes calculations, but it is important to understand their individual effects. 
These will only be briefly discussed here; a more detailed analysis is planned for the future. 

It should be understood that the search for the optimum implicit method should not 
be limited to, and in fact, may very well lead outside the small set of splittings considered 
here. We think of this study only as a pathfinder for this search. 


2. ANALYSIS 


The stationary system of conservation laws for inviscid flows is written as 

OF dG dH 0 „ 

R(u) = & + a^ + ~dl + s - 0 ’ 


(i) 


where U is the vector of conserved state variables, i.e. U = p(l,u,v,w, E) T ] F,G,H are 
the corresponding inviscid fluxes and S is a source term. To make the nonlinear system 
numerically solvable, a Newton linearization procedure is employed, yielding the iterative 
equations 

L(U)6U = — R(U), (2a) 

where 

6U = U'* +1 - U". (26) 

L(U) and R(U) are called the implicit and the explicit operators respectively; the implicit 
operator is given by 


- . dF dG dH \ dR 
L ( U ) ~ dv[ dx + dy + dz + ) ~ dU’ 


( 3 ) 


Spatial discretizations of L(U) and R(U) are made using the upwind-differencing tech- 
niques of the FVS[5,6] and the FDS[7]. 
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For an illustration, let us consider the x-flux F; the other two fluxes follow similarly. 
The Jacobian matrix A of F, 

A = J§, (4.) 

is diagonalizable, with its right eigenvectors forming the columns of the transformation 
matrix Q: 

A = Q _1 AQ, A = diagonal matrix. (46) 

A splitting of A into nonnegative and nonpositive parts yields a splitting of A, 

A* =Q-A ± Q, 

A = A + + A" . 

Steger and Warming[5] constructed the flux splitting 

F = F + + F" 


(5a) 

(56) 

(6a) 


by applying the eigenvalue splitting (5) and exploiting the fact that F is a homogeneous 
function of degree one. Hence the split fluxes and eigenvalues are connected by 

F* = (Q" 1 A ± q)u. (66) 

The associated Jacobian matrices, often called true Jacobians, are denoted by 

dF ± 

(6c) 

Unfortunately this splitting leads to discontinuous switching between F + and F” , that 
manifests itself in “glitches” at sonic and stagnation points, and broadening of shock 
profiles. 

Van Leer [6] sought to remove this flaw by taking an entirely different approach. He 
devised a differentiable splitting of F in the form (6a) by writing the mass flux as the sum 
of two quadratic polynomials in Mach number, and building the momentum and energy 
fluxes on this basis. No requirement of homogeneity of the fluxes is needed. Recently a 
more general derivation has been presented! 11], including a whole family of differentiable 
split fluxes and extending to a general equation of state for real gases. 

Roe[7] constructed the flux-difference splitting by accounting for the wave interaction 
of neighboring cells via an approximate Riemann solution. He connected the difference in 
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fluxes with the difference in state variables by a “mean-value” Jacobian between neighbor- 
ing states, 


-F y _ 1 =A y . 1 /a (U y -U y _ l j, 

(7a) 

A y - 1 /a =A(u y _»,U y ) 


= a(u(u,_ 1 ,u,)) 

(76) 

III 

> 

"S' 

1 

w 



where A is identical in form to A in (4a) but evaluated at the Roe-averaged state U[7]. 
Again, eigenvalue splitting yields matrix splitting, 


A* = Q" l A±Q. 


( 8 ) 


The eigenvalues of A vanish with the characteristic speeds in the x-direction or, in gen- 
eral, in the direction normal to the cell interface. This leads to a crisp representation of 
stationary discontinuity surfaces (shock, contact surface) when these are aligned with the 
interface. 

With the above splittings, the implicit operator L contains the following nonzero 
contributions from the x-flux F: 


FDS : 



The algebra involved in deriving the true Jacobians A* in the FVS is tedious; yet the 
algebra of differentiation in the FDS is formidable. Mulder and Van Leer[l5] simplified 
the algebra considerably by ignoring the spatial derivatives of A ± . The approximate form 
(referred to as ARS) yields the following elements in L: 


ARS 


A + 

- 1/2 5 


(^y-1/2 -^y+1/2)* 1 


/2 ' 


(9c) 


This approximation, as will be seen later, proves to be detrimental in some calculations. 

For the discretization of the residual R let us denote the forward and backward dif- 
ference operators by 

A + (»); = Hy+1 - (•), , 
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A" (•), = (•)> - (•), -!. 

Then the appraxiamtion of dF / dx is given in the compact form: 


3F 


3 _ 


dx 


Ax 


e 


a-f; +a + f; +-a + a- (F;_ 4 -F; +f ) 


+ o 


Ai 1+# ,Ai 2 


(10a) 


(106) 


which includes the following cases: 

6 = 0, first-order 

6 = 1, second-order 

l = 1, full upwinding 

l = 0, central differencing. 

Note that for the evaluation of L(U) in (3) we always assume 0 = 0, while the R(U) 
has the option of taking 0 = 0 or 0 = 1. To eliminate oscillations and obtain a sharp 
representation of discontinuities, we use an upwind-based TVD scheme[l6|. Equation (10b) 
may be modified by introducing a limiter <p multiplying the second-order flux differences. 


WjL 

dx 


1 

Ax 




+ A + FT + *-A+ (<£>y- i A' 






+ 0\Ax 1+e ,Ax 2 

( 10c) 


The limiters are nonlinear functions of neighboring flux differences. For reasons mentioned 
in Section 4 we did not include this in the present study. 


After substitution of spatial discretizations, along with proper treatment of bound- 
ary conditions given in next section, a block matrix system is obtained. To solve this 
system, direct LU decomposition is used for one-dimensional problems and line-relaxation 
procedures are used in the two-dimensional problems. 


3. SPECIAL NUMERICAL PROCEDURES 

In addition to what has been said in general about the numerical methods used, 
two special procedures must be discussed; these regard the entropy condition and various 
boundary conditions. 

§3.1. Satisfying the Entropy Condition 

Of all splitting methods described above, Roe’s is the only one that may violate the 
entropy condition. This happens in particular in the case of a stationary expansion shock 
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located exactly at an interface. Just as for an admissible shock, the Roe matrix A has a 
vanishing eigenvalue. From the general flux formula, 

Fy-i/a — Fy_i + A;., /a (U, - U y . ,) = F, - AJ_ .^(U, - U,_.) 

= i(F,- 1 +F,)-i|A,_ 1 /! |(U,-U,. 1 ). |A| = A + — A" , (U) 

it follows that the flux across the shock equals pre- and post-shock fluxes, and therefore 
can not cause any spreading of the wave. In practice, e.g. in the one-dimensional nozzle 
flows of Section 4.1, this leads to nonuniqueness of the numerical solution in the form of 
an expansion shock located at a sonic point. Its strength is a function of the precise path 
followed in reaching a steady state, and therefore depends on the implicit operator used 
(see Section 1, remark (1)). To overcome this problem one must make sure there is some 
residual dissipation. Following Harten[l7] we do this by modifying the eigenvalues |A| of 
|A| when these are close to zero: 


I A ly- 1/2 * 77 > l A y— 1 / 2 1 < 6 /-i/2> 

4C /-l/2 


with 


Cy_ !/ 2 = K max (Ay — Ay_ l , 0), K > 0. 

This modification is needed only for the genuinely nonlinear characteristic families. 


(12a) 

(126) 


§3.2 Boundary Treatment 

In order to achieve rapid convergence to the steady state, it is necessary to make the 
boundary procedure implicit. Explicit procedures will deteriorate the convergence rate and 
sometimes lead to instability. The question remains whether there is a unique and most 
effective implicit boundary procedure. The answer seems to be that there are perhaps as 
many boundary shemes in the literature as there are interior schemes. It also seems that 
even less definitive appraisals can be made of boundary schemes than of interior schemes. 
This may have the following causes: (l) Boundaries are only a small part of the whole 
domain considered, therefore one paid less attention; (2) There is a strong interaction 
between boundary and interior schemes, both on the level of the discretization and of the 
PDE itself. As Moretti[l8] put it: “Flow evolution depends not only on the events at one 
boundary but on the interaction of the two boundaries and of the interior as well.” The 
interaction is hard to analyze. Often a linear analysis can at best give clues as to the 
stability, but not at all to the convergence rate of the full method. 

To compare the effect of the operators related to different split- flux formulas, it would 
seem reasonable to combine all with the same boundary treatment. This procedure obvi- 
ously conflicts the statement (2) above and therefore may not lead to optimum results for 
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any of the interior schemes considered here. However, this is the best we can do, judging 
the present state of analysis. To illustrate our implicit boundary procedure, we consider 
subsonic outflow from a boundary point, denoted by subscript J. The static pressure is 
fixed, i.e., 

Pj = const, (13a) 

or 

6pj = 0. (136) 

This reduces to 

S(pE)j — (uSpu - -u?6p) . (13c) 

Allowing extrapolation of two variables from interior points, denoted by J — 1 and J — 2, 
we choose 



(a + 1) 


Sp \ 

SpvJj-i 


a = 0 or 1. 


For the equation centered 


at J — 1, the operator L(U) has a term 


Aj SVj 


( Oil 
021 
o 3 i 


O 12 

O 22 

O 32 


Ol3 \ 

023 I 


a 33 J j 



( O11X o 12 /? 
O 21 X a 22^ 
O31X O32 /? 

( OllX 
O21X 
a 3 i X 

( OnX 
O21X 
O31X 



di 2 ft o\ ( 6 p \ 

a ?20 0 I j 6 pu I 

o 32 /3 0 J j \ 6 pE J J _ l 

012/5 0 \ / Sp \ 

022/? 0 I 6 pu j 

0-32 0 0/ J J J_ 2 


(14) 


(15) 


where x = (l — 2 u2 ) an ^ = 1 + u; the elements (au,...,a 33 ) of A - are known at J 
at the last iteration level. By substitution of (15) in the discrete equation centered at 
J — 1, the boundary term Aj 6U j is now expressed in the interior terms |A./_ x |6U./_ ! 
and Aj_ 2 6U/_ 2 - A similar procedure can be performed for a subsonic inflow boundary 
where we specify, as suggested by Yee et al[19], 

Pi — const, 
u x = const. 


(16a) 
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Hence a spatial extrapolation to interior points gives 

6pE j - l) _1 [(a + 1 )Sp 2 - a^p 3 ] 

= {a + l)[6pE -uSpu + ^Bp}^ (166) 

— a\6pE — uSpu + ^u 2 £p] 3 , 

£t 

and Af 5U! yields contributions to the terms associated with £U 2 and 6U 3 . 

As for supersonic inflow one simply sets 

ev, = o, (n) 

and at the supersonic outflow boundary, 

6U j = (or + l) 5 U/_ i — 06U j _ 2 • ( 18 ) 


4. NUMERICAL TESTS 


Numerical tests of the various combinations of implicit and explicit operators have 
been carried out for the Euler equations in one and two dimensions, and for the Navier- 
Stokes equations in two dimensions. The one-dimensional tests are nozzle-flow problems 
that are governed by the equations 

— - + S = 0, (19a) 

OX 

where the source term S includes the effect of nozzle area s(x) 



$(x) = s'(x)/s(x). 


The area distribution for the divergent nozzle is given by 

s(x) = 1.398 + 0.347 tanh( 0.8z - 4), 0 < x < 10, 


and for the convergent-divergent nozzle by 

s(x) = 1.75 — 0.75 cos(x — 5)7r/5, 0 < z < 5 
s(x) = 1.25 — 0.25 cos(x — 5)n/5, 5 < z < 10. 


(196) 


(20a) 


(206) 
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Unless noted, 100 uniform meshes were used in the one-dimensional calculations. 

The two-dimensional Euler problem concerns the reflection of a Mach 2.9 oblique 
shock wave off a solid surface at y — 0; the angle /? between incident wave and wall is 29 
degrees. The computational domain, 0 < x < 4 and 0 < y < 1, was equally divided into 
60 by 20 meshes. 

Finally, the Navier-Stokes calculation concerns the interaction of an oblique shock 
wave with a laminar boundary layer; the flow parameters are those of recent experi- 
ments^©] and are given later, along with details about the computational grid used. 

§4.1 One-Dimensional Problems 

For the one-dimensional nozzle problems, the variations in combinations of operators 
lead only to minor difference in accuracy and convergence rate if the flow is either fully 
supersonic or subsonic; these results therefore are not discussed in detail. Essential differ- 
ences appear as a shock wave enters in the solution; the results are summarized in Tables 1 
to 5 and Figures 1 to 4. Except for the changes in the flux formulas all calculations of the 
same flow were done for the same set of numerical parameters ( grid spacing, relaxation 
factor) and with the same boundary treatment.This does not imply that any combination 
of operators achieves its peak performance. At this stage, our goal is merely to select a 
pair of operators that maintains stability, accuracy and reasonably fast convergence (say, 
within 100 iterations) for a wide range of problems. Once this goal has been reached one 
may start looking for ways to further improve the efficiency of the method. 

Convergence is measured by summing up all changes in the state quantities: 

< 21 ) 

grids t= 1 

for a solution to be called “converged” e must have been reduced to 0(10“ 11 ). The first set 
of results is for divergent-nozzle flow. Table 1 shows the error and number of iterations till 
convergence for all possible first-order and second-order accurate upwind approximation of 
R(U); for L(U) only the first-order approximations are used. To assure that convergence 
will be achieved at least for all first-order methods, an under-relaxation factor — varying in 
the range of 0.15 to 0.6 — was applied in all Newton iterations. It is seen that the fastest 
convergence rate is generally achieved using Van Leer’s splitting in L(U), with approximate 
Roe’s splitting (ARS) coming in second for first-order calculations of R(U) but becoming 
totally useless in second-order calculations. In the latter case the iterations are led into a 
limit cycle. 

The number of iterations needed appears to be lowest, regardless of the choice splitting 
in L(U), when R(U) is discretized with SWS. This may have to do with the fact that SWS 
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is the most dissipative of the three splittings considered. This dissipative effect is evident 
from Fig. 1; the SWS shock structure is clearly smeared on the subsonic side. The RS 
and VLS results are both excellent. Note that the results on the supersonic side of the 
shock are identical for all schemes, since the schemes are identical in that region and the 
solution is completely determined by the inflow conditions. Quantitative results are found 
in Table 2, where numerical values of the pressure, normalized by the point values of the 
exact solution, are listed. 

In the second-order solutions of Fig.2 all flux formulas have produced a large overshoot 
of the Mach number in the pre-shock state; for SWS there is also a tendency to undershoot 
the post-shock state. This non-monotonicity could have been prevented by the introduction 
of limiters in the flux-differencing formula, as in (10c) (see [16] for details). It was found 
though that these limiters — nonlinear by design — have a strong effect on the convergence, 
always slowing it down and sometimes ending it in a limit cycle, although the error at 
which the limit cycle happens is far below the local discretization error. It thus appears 
that the effect of limiters should be studied separately; they have not been used in any 
of the one-dimensional calculations. Regardless of the performance near the shock, all 
second-order discretizations give better accuracy away from the shock than the first-order 
discretization, as evident from Table 2. 

The convergent-divergent nozzle used for the next series of tests has a discontinuity in 
curvature at the throat, leading to a discontinuity in the first derivative of state variables. 
This makes it a good test of the accuracy of splittings. From Fig. 3 and 4 it is evident 
that in the area of accuracy, i.e. in R(U), RS and VLS perform much better than SWS, 
which introduces an O(Ax) jump at the sonic point. 

Table 3 summarizes the convergence histories and indicates the need for consistency 
of L(U) and R(U). Invariably, a degradation of efficiency is observed when splittings 
are mixed. The worst matches are seen for SWS in R(U), which can only be made to 
converge using SWS also in L(U). In contrast, VLS and ARS/RS appear to be more or 
less interchangeable. While far from optimal, especially in the second-order case, SWS 
still seems to be most generally applicable splitting for use in L(U). On the other hand, if 
one would preclude SWS in R(U) for its inaccuracy, the choice for L(U) clearly belongs 
to VLS. 

In order to further study the applicability of SWS as a multi-purpose splitting in 
L(U) we removed the under-relaxation factor for the Newton iteration, allowing each 
scheme to converge at the maximum possible rate. With the same splitting and the same 
accuracy in L(U) and R(U), quadratic convergence may be achieved. The results for 
the divergent nozzle are listed in Table 4. Quadratic convergence is indeed observed for 
SWS/SWS and VLS/VLS with first-order accuracy; these are the only combinations that 
form a true Newton method. The general applicability of SWS in L(U) is evident; without 
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under-relaxation VLS or ARS in L(U) can not be mixed with another splitting in R(U). 
Second-order accuracy with RS in R(U) can not even be obtained with ARS as the implicit 
driver, while SWS accomplishes this without any problem(see remark (2) in Section 1). 

A final experiment, this time for the convergent-divergent nozzle, was to double Ax 
and observe the effect on convergence and stability. The results are compiled in Table 5; 
these show no significant changes with respect to those in Table 3. 

From these one-dimensional experiments we may draw some provisional conclusions. 

Regarding the implicit operator: 

(1) SWS is prefered over VLS and ARS as a general technique, because of its robustness. 
It is applicable, even without under-relaxation, with each of the splittings used in the 
explicit operator, and achieves convergence — even with second-order approximation 
of R(U) — within 100 iterations. 

(2) Implicit VLS combined with explicit VLS leads to the best convergence rate, even for 
the second-order explicit operator. 

(3) ARS in the form (9c) is not an acceptable implicit discretization, not even in com- 
bination with explicit RS. The subject of approximating the Jacobian of Roe’s flux 
clearly needs futher study. 

Regarding the explicit operator: 

(1) RS and VLS are far superior to SWS regarding the numerical accuracy achieved; it 
may be noted that VLS requires fewer operations than RS. 

From the above conclusions one may further derive that VLS used in both implicit 
and explicit operators is a “best buy” as to the accuracy and convergence rate achieved in 
the one-dimensional nozzle problems. When comparing the above conclusions with (i) and 
(ii) in Section 1, we have not yet been able to make a clear distinction in accuracy between 
VLS and RS. This is due to the absence of a contact discontinuity in the one-dimensional 
test problems. This issue will be resolved by the Navier-Stokes experiment of Section 4.3. 

§4.2 Two-Dimensional Euler Problem 

For the shock-reflection problem we may allow ourselves to concentrate on accuracy. 
Owing to the flow being fully supersonic in the x-direction, convergence is achieved within 
eight iterations for all schemes using line/Gauss-Seidel relaxation with the line in the y- 
direction and downstream sweeping in the x-direction. The other side of the coin is that 
the upwind schemes lose their property of representing a shock front by a narrow profile. 
This is illustrated by Fig. 5, showing the pressure distribution on the wall (y = 0) and at 
the mid-plane (y = 0.5), with RS and first-order accuracy. When going to second-order 
accuracy, resolution is greatly improved, as shown in Fig. 6; there is hardly any difference 
between using RS or VLS in R(U). (SWS was dropped as a discretization for R(U) 



because of its inferior performance in the tests of Section 4.1) The pressure undershoots 
and overshoots near the shocks can be avoided by using the limiters, see (10c), as is 
illustrated by Fig. 7. There is a large choice of limiters[2l]; here we used the “minmod” 
function for illustration. Wall pressures from the first-order and second-order solutions are 
also tabulated in Table 5. 

Although other limiters — for example “superbee” — may sharpen the shock profile, it 
should be noted that any genuine improvement of the representation of oblique shocks 
in supersonic flow must involve either grid adaptation or some kind of wave-recognition 
algorithm [22]. 

As anticipated, the results of the above two-dimensional Euler tests do not add any- 
thing to the results of the one-dimensional tests in regard to discriminating between VLS 
and RS. 


§4.3 Two-Dimensional Navier-Stokes Problem 


As for the Navier-Stokes calculations, we again concentrate on the issue of accuracy. 
The SWS is assumed an acceptable driver in L(U) while R(U) is varied according to the 
three flux formulas. The thin-layer form of the Navier-Stokes equations results in these 
operators: 


% dF dG dW 
R < u ' u »> = ^ + ^--^ = 0 ' 


( 22 a) 


where F and G are again the inviscid fluxes and W is the viscous flux that involves only 
the y-derivatives 


W(U,U„) r 


7 M 2 M L du 4 dv du 

1 °° fz 0 U h 

Re oo ’ dy ’ 3 dy ’ dy 


Av dv 7 de 
3 dy + P r dy 


and 


L < u > = ^llj ro 




a 2 

a y 2 



(226) 


A detailed description of the derivation and the resulting algebraic system can be found 
in [23], While the convective terms continue to be upwind-differenced, but the diffusion 
terms are now central-differenced because the latter have no bias in direction. 


To minimize the uncertainties regarding turbulence modeling in the tests, we restrict 
ourselves to laminar flow. The problem chosen for tests is the interaction of an oblique 
shock wave with a boundary layer over a flat plate, recently reported in [20]. The flow 
parameters are: = 2.2, R eoo — l.Ox 10 5 , and the wedge angle of 3.75 degrees producing 

a shock angle of 30.027 degrees; the reference length used in Re «, is the distance from the 
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leading edge to the point of shock impingement on the plate. The wall is adiabatic, the 
Prandtl number is 0.72 and the viscosity is assumed to follow Sutherland’s formula. 

It was found in [23] that a grid distribution deduced from the triple-deck theory[24] 
gave results more accurate than those not based on such a grid [25]. Here we again use the 
same spacing, 

r Ax = 0.03, 

Ay, = 1.5625e(-4), 1 < t < 4, 

Ay, = 1.1868Ay,_i, 5 < * < 33, 

, Ay, = 3.75e(— 2), i > 33; 

the total number of meshes is 74x64. The co-ordinates x and y have been non-dimensionalized 
by the same reference length for Re*, • 

The flow remains supersonic on all boundaries except at and near the wall. The 
computational domain is chosen such that the reflected shock wave from the wall does 
not intersect the top boundary and the flow is sufficiently well-developed after leaving the 
interaction region. Hence it is appropriate to specify all variables at the inflow and the 
top boundaries. Linear extrapolation from the interior points is used for the variables at 
the outflow. 

In [23] only SWS was used; at present we also include the other splittings for com- 
parison. Figure 8 shows the calculated and measured surface pressure; Roe’s splitting 
is clearly more accurate than the Van Leer and Steger- Warming splittings which tend 
to diffuse details of the solution (see Section 1, remark (4)). Convergence also happens 
to be fastest with RS, taking about 250 iterations as compared to 400 for VLS and 520 
for SWS for e to be reduced by five orders of magnitude, although Roe’s splitting formula 
needs slightly more computations in each iteration. It is worth mentioning that the present 
upwind-differencing results are much closer to the data than the central-differencing results 
in [20]. 


5. CONCLUSIONS 

The present study confirms some previously known results regarding the accuracy of 
various splitting algorithms, and adds some regarding their use in implicit operators. As to 
the question of accuracy, Roe’s flux-difference splitting clearly is superior to the two flux- 
vector splittings considered — Van Leer’s and Steger- Warming’s — and should be considered 
the best choice for use in a Navier-Stokes code. If only Euler calculations are considered 
and contact discontinuities or vortex sheets are not important, Van Leer’s splitting is an 
attractive alternative. All splittings give insufficient resolution of an oblique shock wave 
in supersonic flow. 
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As to convergence, the Steger- Warming splitting is the only splitting of the implicit 
operator that can be combined with any of the splittings of the explicit operators and may 
therefore called the most robust of all splittings. This conclusion is of limited importance, 
since accuracy requirements clearly exclude the use of the Steger- Warming splitting in 
the explicit operator. Under this constraint, the Van Leer splitting becomes the better 
splitting for the implicit operator, since it leads to faster convergence with the remaining 
splittings of the explicit operator. The data, however, are inconclusive, since a Navier- 
Stokes computation with the Van Leer splitting as the implicit operator is yet missing. 

Least effective on the implicit side is the approximate Roe splitting, which in second- 
order calculations even fails to lead to convergence in combination with the Roe splitting 
on the explicit side. It must be mentioned that the approximate Roe splitting, precisely 
because it is approximate, is not unique; many variations may have to be tested sepa- 
rately.(Barth[13] considers two and explains their different behavior in earlier and later 
stages of convergence.) This situation is far from desirable, but the use of the exact Ja- 
cobian of the Roe splitting in the implicit operator is even less attractive, because of its 
great computational cost. 

Future extensions of the present study will concentrate on the question of selecting for 
the discretized Navier-Stokes equations a splitting of the implicit operator that is efficient 
and robust when combined with the Roe splitting of the explicit (residual) operator. 
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Figure 1. Divergent nozzle, Woo — 1-26, Pe/Pt = 0.746; 
first-order accurate R(TJ). Solid line: exact 
solution, symbols: numerical solution. 


* 



Figure 2. Divergent nozzle, Moo = 1-26, P e /Pt = 0.746; 

second-order accurate R(U) . Solid line: exact 
solution, symbols: numerical solution. 


Table 1 

Convergence of numerical solutions for a divergent nozzle flow; 
Woo = 1.26 y P e /P t - 0.746. The maximum number of 
iteration allowed was 300. 


L(U) 

R(U) 

O(Az) T 

0(Az 2 ) ~ 

Error 

No. of iter. 

Error 

No. of iter. 

sws 

SWS 

6.692(-ll) 

67 

6.635(-ll) 

112 

VLS 

7.828(-ll) 

118 

8.230(-ll) 

127 

RS 

8.786(-ll) 

120 

8.576(-ll) 

127 

VLS 

SWS 

6.676(-ll) 

68 

6.297 (- li) 

106 

VLS 

8.836(-ll) 

75 

7.233(-ll) 

105 

RS 

7.929(-ll) 

81 

8.881(-11) 

105 

ARS 

SWS 

7.262 (-11)^ 

75 


— 

VLS 

9.087(-ll) 

93 

— 

— 

RS 

7.325(-ll) 

94 

— 

— 


— Stagnant or slow convergence 
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Table 2 

Ratio of numerical to exact solutions for a divergent nozzle flow; 
Moo = 1.26, P e /P t =0.746. 



Pcomp/ 

Pczact 

X 

0( Ax) 

O(Ax^) | 


sws 

VLS 

RS 

SWS 

VLS 

RS 

4.0 

4.1 

4.2 

4.3 

4.4 

4.5 

4.6 

4.7 

4.8 

4.9 

5.0 

5.1 

5.2 

5.3 

5.4 

5.5 

5.6 

5.7 

5.8 

5.9 
6.0 

0.98583(0) 

0.98528(0) 

0.98489(0) 

0.98467(0) 

0.98465(0) 

0.98485(0) 

0.98527(0) 

0.10445(1) 

0.21495(1) 

0.80460(0) 

0.88994(0) 

0.93628(0) 

0.95951(0) 

0.97113(0) 

0.97740(0) 

0.98126(0) 

0.98401(0) 

0.98620(0) 

0.98803(0) 

0.98962(0) 

0.99101(0) 

0.98583(0) 

0.98528(0) 

0.98489(0) 

0.98467(0) 

0.98465(0) 

0.98485(0) 

0.98527(0) 

0.98593(0) 

0.15310(1) 

0.94106(0) 

0.98294(0) 

0.98412(0) 

0.98534(0) 

0.98657(0) 

0.98779(0) 

0.98897(0) 

0.99011(0) 

0.99118(0) 

0.99218(0) 

0.99310(0) 

0.99394(0) 

0.98583(0) 

0.98528(0) 

0.98489(0) 

0.98467(0) 

0.98465(0) 

.0.98485(0) 

0.98527(0) 

0.98593(0) 

0.14962(1) 

0.98479(0) 

0.98751(0) 

0.98850(0) 

0.98947(0) 

0.99043(0) 

0.99135(0) 

0.99223(0) 

0.99307(0) 

0.99385(0) 

0.99456(0) 

0.99522(0) 

0.99581(0) 

0.99944(0) 

0.99953(0) 

0.99964(0) 

0.99977(0) 

0.99992(0) 

0.10001(1) 

0.10003(1) 

0.66515(0) 

0.17969(1) 

0.94302(0) 

0.10384(1) 

0.10241(1) 

0.10035(1) 

0.99769(0) 

0.99851(0) 

0.99965(0)' 

0.99995(0) 

0.99991(0) 

0.99987(0) 

0.99987(0) 

0.99988(0) 

0.99944(0) 

0.99953(0) 

0.99964(0) 

0.99977(0) 

0.99992(0) 

0.10001(1) 

0.10003(1) 

0.81497(0) 

0.14880(1) 

0.10412(1) 

0.10144(1) 

0.99985(0) 

0.99980(0) 

0.99979(0) 

0.99978(0) 

0.99979(0) 

0.99980(0) 

0.99981(0) 

0.99983(0) 

0.99984(0) 

0.99986(0) 

0.99944(0) 

0.99953(0) 

0.99964(0) 

0.99977(0) 

0.99992(0) 

0.10001(1) 

0.10003(1) 

0.82237(0) 

0.14924(1) 

0.10530(1) 

0.10002(1) 

0.99979(0) 

0.99982(0) 

0.99981(0) 

0.99980(0) 

0.99980(0) 

0.99981(0) 

0.99982(0) 

0.99983(0) 

0.99984(0) 

0.99986(0) 
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Figure 3. Convergent-divergent nozzle, Moo = 0,2395, 

P e /P t = 0.84; first-order accurate R(U). Solid line: 
exact solution, symbols: numerical solution. 



Figure 4. Convergent-divergent nozzle, Moo = 0.2395, 

Pe/Pt = 0.84; second-order accurate R(U). Solid 
line: exact solution, symbols: numerical solution. 


Table 3 

Convergence of numerical solutions for a convergent-divergent 
nozzle flow; M^ = 0.2395, P e /P t = 0.84. 


L(U) 

R(U) 

0( Ax) 

0(Ax 2 ) 

Error 

No. of iter. 

Error 

No. of iter. 

sws 

SWS 

8.4191-11) 

189 

8.325(-ll) 

187 

VLS 

2.068 (-06J 

300 

1.695 (-06) 

300 

RS 

3.030(-10) 

300 

1.944(-06) 

300 

VLS 

SWS 

— 

— 

— 

— 

VLS 

8.542(-ll) 

204 

9.034(-ll) 

176 

RS 

6.857(-08) 

300 

7.002(-ll) 

175 

ARS 

SWS 

— 

— 

— 

— 

VLS 

4.661f-08) 

300 

6.9221-08) 

300 

RS 

8.719(-ll) 

189 

9.290(-ll) 

200 


— Stagnant or slow convergence 
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Table 4 

Improvement of convergence rate for a divergent nozzle flow; 
Moo = 1.26, P e /P t = 0.746. 


L(U) 

R(U) 

O(Ax) 

0(Az 2 ) 

Error 

No. of iter. 

Error 

No. of iter. 

SWS 

SWS 

8.464(-ll) 

23 

8.911(-11) 

63 

VLS 

7.688 (-11) 

45 

5.139(-ll) 

79 

RS 

6.907(-ll) 

45 

6.776(-ll) 

79 

VLS 

VLS 

5.594(-ll) 

23 

5:682 (-11) 

59 

ARS 

RS 

6.846(-ll) 

46 

— 

- - 


Unstable, also for other pairs not listed. 


Table 5 

Effect of doubling Ax on convergence of numerical solutions for 
a convergent-devergent nozzle flow; M a 0 = 0.2395, P e /P t = 0.84. 


L(U) 

R(U) 

0( Ax) 

0(Ax 2 ) 



Error 

No. of iter. 

Error 

No. of iter. 


SWS 

9.118f-ll) 

159 

9.280(-ll) 

211 

SWS 

VLS 

4.012(-06) 

300 

3.082(-06) 

300 


RS 

1.987 (-06) 

300 

1.180(-06) 

300 


SWS 

— 

— 

— 

— 

VLS 

VLS 

7.324(-ll) 

192 

rH 

1 

s 

00 

od 

176 


RS 

8.289(-l 1) 

189 

9.843 (-11) 

179 


SWS 

— 

— 

— 

— 

ARS 

VLS 

1.962 (-07) 

300 

3.76lf-07] 

300 


RS 

1.601 (-06) 

300 

1.180(-11) 

221 


Stagnant or slow convergence 
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Figure 5. Shock reflection problem, = 2.9, fi — 29°; 

first-order accurate R(TJ) using R S; I ++: exact 
solution, o o o: numerical solution. 
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Figure 6. Shock reflection problem, = 2.9, = 29°; 

second-order accurate R(U) using (a) VLS and (b) 
RS; -H H-: exact solution, ooo: numerical solution. 





p/i 








Table 6 

Comparison of numerical and exact solutions for a 
shock reflection problem; Moo = 2.90,^ = 29.0°. 


X 

vail 

Exact 

0[ Ax 2 ) 

TVD 

VLS | RS 

VLS RS 

1.000 

1.067 

1.133 
1.200 

1.267 

1.333 

1.400 

1.467 

1.533 
1.600 

1.667 

1.733 
1.800 

1.867 

1.933 
2.000 

2.067 

2.133 
2.200 

2.267 

2.333 

2.400 

2.467 

2.533 
2.600 

2.667 

2.733 
2.800 

2.867 

2.933 
3.000 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.4108(1) 

0.1007(1) 

0.1012(1) 

0.1013(1) 

0.1006(1) 

0.9843(0) 

0.9447(0) 

0.8886(0) 

0.8578(0) 

0.9534(0) 

0.1252(1) 

0.1742(1) 

0.2342(1) 

0.2936(1) 

0.3429(1) 

0.3778(1) 

0.3988(1) 

0.4095(1) 

0.4138(1) 

0.4147(1) 

0.4144(1) 

0.4137(1) 

0.4130(1) 

0.4125(1) 

0.4122(1) 

0.4121(1) 

0.4119(1) 

0.4118(1) 

0.4116(1) 

0.4114(1) 

0.4111(1) 

0.4109(1) 

0.1007(1) 

0.1012(1) 

0.1015(1) 

0.1011(1) 

0.9912(0) 

0.9501(0) 

0.8862(0) 

0.8430(0) 

0.9343(0) 

0.1237(1) 

0.1725(1) 

0.2309(1) 

0.2885(1) 

0.3370(1) 

0.3720(1) 

0.3938(1) 

0.4052(1) 

0.4099(1) 

0.4113(1) 

0.4114(1) 

0.4112(1) 

0.4112(1) 

0.4114(1) 

0.4116(1) 

0.4118(1) 

0.4120(1) 

0.4120(1) 

0.4119(1) 

0.4117(1) 

0.4114(1) 

0.4111(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1001(1) 

0.1003(1) 

0.1010(1) 

0.1048(1) 

0.1248(1) 

0.1744(1) 

0.2353(1) 

0.2896(1) 

0.3425(1) 

0.3766(1) 

0.4005(1) 

0.4065(1) 

0.4094(1) 

0.4110(1) 

0.4117(1) 

0.4120(1) 

0.4120(1) 

0.4119(1) 

0.4118(1) 

0.4117(1) 

0.4115(1) 

0.4113(1) 

0.4111(1) 

0.4108(1) 

0.4105(1) 

0.4102(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1000(1) 

0.1001(1) 

0.1004(1) 

0.1031(1) 

0.1246(1) 

0.1750(1) 

0.2324(1) 

0.2852(1) 

0.3392(1) 

0.3761(1) 

0.4045(1) 

0.4070(1) 

0.4077(1) 

0.4090(1) 

0.4103(1) 

0.4106(1) 

0.4110(1) 

0.4111(1) 

0.4110(1) 

0.4108(1) 

0.4104(1) 

0.4103(1) 

0.4103(1) 

0.4104(1) 

0.4104(1) 

0.4104(1) 
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P/Poo P/A 




Figure 8. Surface pressure distribution P/Poo for a shock- 
wave /laminar boundary- layer interaction, 

A4o = 2.2, ikoo = 1.0 x 10 5 and 0 = 30.027°; 
second-order accurate R(U) using (aj RS, (b) VLS 
and (c) SWS; symbols: measurement [20], solid line: 
computation. 
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